Crystalline Particle Packings on a Sphere with Long Range Power Law Potentials 
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The original Thomson problem of "spherical crystallography" seeks the ground state of electron 
shells interacting via the Coulomb potential; however one can also study crystalline ground states 
of particles interacting with other potentials. We focus here on long range power law interactions 
of the form 1/r 7 (0 < 7 < 2), with the classic Thomson problem given by 7 = 1. At large R/a, 
where R is the sphere radius and a is the particle spacing, the problem can be reformulated as 
a continuum elastic model that depends on the Young's modulus of particles packed in the plane 
and the universal (independent of the pair potential) geometrical interactions between disclination 
defects. The energy of the continuum model can be expressed as an expansion in powers of the 
total number of particles, M ~ (R/a) 2 , with coefficients explicitly related to both geometric and 
potential-dependent terms. For icosahedral configurations of twelve 5-fold disclinations, the first 
non-trivial coefficient of the expansion agrees with explicit numerical evaluation for discrete particle 
arrangements to 4 significant digits; the discrepancy in the 5th digit arises from a contribution to 
the energy that is sensitive to the particular icosadeltahedral configuration and that is neglected in 
the continuum calculation. In the limit of a very large number of particles, an instability toward 
grain boundaries can be understood in terms of a "Debye-Huckel" solution, where dislocations have 
continuous Burgers' vector "charges". Discrete dislocations in grain boundaries for intermediate 
particle numbers are discussed as well. 

PACS numbers: 61.72.Mm,61.72.Bb,64.60.Cn,82.70.Dd 



I. INTRODUCTION 



The Thomson problem of constructing the ground 
state of (classical) electrons interacting with a repulsive 
Coulomb potential on a sphere Q has a rich, approxi- 
mately one hundred year old history 0- S E3| ■ Determin- 
ing crystalline particle packings in curved geometries has 
a number of interesting applications in physics, mathe- 
matics, chemistry and biology particularly if one allows 
more general interactions amongst the particles. 

An almost literal realization of the Thomson problem 
is provided by multi-electron bubbles 0, @ . Electrons 
trapped on the surface of liquid helium by a submerged, 
positively charged capacitor plate have longbeen used 
to investigate two dimensional melting 0, 13- Multi- 
electron bubbles result when a large number of electrons 
(10 5 — 10 7 ) at the helium interface subduct in response to 
an increase in the anode potential and coat the inside 
wall of a helium vapor sphere of radius 10—100 microns. 
Typical electron spacings, both at the interface and on 
the sphere, are of order 2000 Angstroms, so the physics is 
entirely classical, in contrast to the quantum problem of 
electron shells which originally motivated J.J. Thomson 
PJ. Information about electron configurations on these 
bubbles can, in principle, be inferred from studying capil- 
lary wave excitations [9|- Similar electron configurations 
should arise on the surface of liquid metal drops confined 
in Paul traps lQj. 

A Thomson-like problem also arises in determining the 
arrangements of the protein subunits which comprise the 



shells of spherical viruses [Til IT^. Here, the "particles" 
are clusters of protein subunits arranged on a shell. Al- 
though the proteins interact predominantly with short 
range Van der Waals potentials, the same issues of spher- 
ical crystallography arise in these protein shells as in the 
original Thomson problem. In spherical viruses, 12 of 
these protein clusters sit at the vertices of a regular icosa- 
hedron in a 5-fold symmetric environment. The remain- 
ing "particles" have 6 neighboring clusters. This problem 
of protein arrangements was solved in a beautiful paper 
by Caspar and Klug f° r intermediate values of R/a, 
where R is the sphere radius and a is the mean parti- 
cle spacing. Caspar and Klug constructed icosadeltahe- 
dral particle packings characterized by integers P and Q, 
which provide regular tessellations of 



M = 10 {P 2 + PQ + Q 2 



(1) 



protein clusters, or "particles" , on the sphere. Most 
known viruses (examples with M as large as 1472 are 
known 0, llil Hlj l fall into this classification scheme, 
and can be studied by use of the continuum methods 
discussed in this paper ^(|. The Caspar-Klug tessella- 
tions of the sphere provide an excellent starting point for 
finding low energy particle configurations on the sphere 
for intermediate values of M w ^-(R/a) 2 . Particles num- 
bers M not in the form of Eq.QJ can be accommodated 
by introducing vacancies or interstitials into these regular 
packings (see |l7| for a discussion of vacancy and inter- 
stitial energies with power law potentials in flat space). 
New ground states involving grain boundaries are needed, 
however, for M > M c w 400 — 600, and in particular in 
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the limit M —> oo [H [U El E3 

Other realizations of Thomson-like problems include 
regular arra ngem ents of colloidal particles in "colloido- 
some" cages |23, H^T, proposed for protection of cells 
or drug-containing vesicle s |25l and fullerene patterns of 
carbon atoms on spheres |2|| and other geometries. An 
example with long range (logarithmic) interactions is pro- 
vided by the Abrikosov lattice of vortices which would 
form at low temperatures in a superconducting metal 
shell with a large monopole at the center [2fl |. In prac- 
tice, the "monopole" could be approximated by the tip 
of a long thin solenoid. 

The problem of best possible packing on spheres has 
also applications in the micropatterning of spherical par- 
ticles [2|| relevant for photonic crystals or Clathrin cages, 
responsible for the vesicular transport of cargo in cells 
[23 (see [3(J for a detailed theoretical study). Crystalline 
domains covering a fraction of the sphere are also of ex- 
perimental interest. In the context of lipid rafts [3l| . 
confocal fluorescence microscopy studies have revealed 
the coexistence of fluid and solid domains on giant unil- 
amellar vesicles made of lipid mixtures. The shapes of 
these solid domains include stripes of different widths 
and orientations The application of spheri- 

cal elasticity to predict shapes of lipid mixtures domains 
has been discussed in 123, |36j • 

In the continuum approach used here, details associ- 
ated with different particle interactions for the system 
discussed above are parameterized by a bulk and shear 
elastic constant and a defect core energy. In practice, 
defect patterns involving dislocations and disclinations 
depend only on the Young's modulus and a core energy 
|2J, which can be determined from flat particle configu- 
rations. Although we concentrate on the computationally 
challenging problem of long range power law potentials, 
explicating and complementing previous results [37|], it 
would be straightforward to treat short range potentials 
as well 0. 

The organization of the paper is as follows. In Sect.lTTl 
some known results for crystals on curved surfaces are 
reviewed and several new results are obtained. The free 
energy of the system is described in Sect. IIIII The par- 
ticular case of the sphere, the Thomson problem, is dis- 
cussed in Sect. IIVI and several predictions for spheri- 
cal lattices with icosahedral symmetry are obtained and 
compared with the results of direct minimizations of dis- 
crete icosadeltahedral particle arrays. The solution of the 
Thomson problem for a very large number of particles is 
discussed in Sect. Sect. IVII contains a summary and 
conclusions, and some technical results are discussed in 
the appendices. 



II. CRYSTALS OF POINT PARTICLES 

Consider a collection of classical point particles con- 
strained to a frozen (non-dynamical) two-dimensional 
surface K, embedded in three-dimensional Euclidean 



space. The particles interact through a general poten- 
tial defined in the three dimensional embedding space or 
solely within the 2d curved surface itself. This paper 
focuses primarily on the potential 



V(R) 



(2) 



Here, e is an "electric charge" such that if R is some 
quantity with dimensions of length, 



e 2 /i? 7 = dimension of energ; 



y 



The case 7 = 1 corresponds to the Coulomb potential in 
three dimensions. Allowing we do not treat this problem 
explicitly here, the replacement 



V(R)^-(\Rr-l) , 

7 



(3) 



allows us to treat the two dimensional Coulomb potential 
by taking the limit 7 — > 0, 



V(R)^-e 2 \og(\R\ 



(4) 



The electrostatic energy of a system of M particles at 
positions interacting via Eq.©, with 1 = (Zi , Z2) , 

Zi, I2 E Z, becomes 



M 

2E = y — — 



R(V) 



(5) 



Note that with this definition the power law interaction 
acts across a cord of the sphere, as would be the case 
for electron bubbles in helium. The focus of this paper 
is the study of crystals on curved surfaces, in particular 
spherical crystals. There are, however, some quantities 
which are insensitive to the curvature of the surface, and 
the simpler geometry of the plane can be used to com- 
pute them. The following two subsections hence focus on 
planar crystalline arrays of particles interacting via the 
potential Eq.©. 



A. Planar Crystals 

The electrostatic energy Eq. (JjJ and the corresponding 
elastic tensor, from which follows the elastic constants of 
the system, may be explicitly computed for crystalline 
orderings of particles in a triangular lattice. 

For any non-compact surface /C, like the plane, the en- 
ergy Eq.JSJ) is divergent for all 7 < 2. If 7 > 0, the 
divergence comes exclusively from the zero mode G = 
associated with the thermodynamic limit of infinite sys- 
tem size. This term (which would be subtracted off if a 
uniform background charge were present) can be isolated 
by setting G = e <C 1 for this mode. The detailed calcu- 
lation is somewhat involved and is given in Appendix lAl 
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The final result for the ground state energy reads 



2£n 



Me 



Me 2 7T x 

2 * r(i 



7(2-7) 



cr( 7 ) 



2j lim 



A C 7/2 GUe |G| 2 ^ 

7 /2 



(6) 
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2 2- 7 
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7/2 gU E IGI 2 -^ 



Ac is the area of the unit cell of the triangular Bravais 
lattice {Ac = ^a 2 ) and V is the Euler Gamma function. 
The coefficient a is a sum over Misra functions, defined in 
Ea. (IB2|) of Appendix IbI The coefficient 0( 7 ) parameter- 
izes the nonsingular part of the energy; its dependence 
on the exponent 7 is shown in Table [I] This negative 
quantity parameterizes the binding energy of the trian- 
gular lattice after the positive "zero mode" contribution 
is subtracted off. For 7 = 1, we have a two dimensional 
"jcllium" model. In the problem considered in the in- 
troduction, no neutralizing background is present, and 
the energy is rendered finite by restricting the crystal to 
a compact surface, like the sphere. The maximum dis- 
tance between points in the surface will then provide an 
infrared cut-off. 

For small displacements of the particles from their 
equilibrium positions, one has 
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1.875 


0.699652 


31 


47.763 


0.875 


0.199772 


23/9 


3.2471 


1.75 


0.619256 


15 


22.647 


0.75 


0.159010 


11/5 


2.7138 


1.625 


0.544152 


87/9 


14.288 


0.625 


0.122622 


21/11 


2.283 


1.5 


0.474268 


7 


10.118 


0.5 


0.090439 


5/3 


1.9294 


1.375 


0.409548 


27/5 


7.625 


0.375 


0.062279 


1.46154* 


1.6352 


1.25 


0.349812 


13/3 


5.9701 


0.25 


0.037955 


9/7 


1.3881 


1.125 


0.295033 


25/7 


4.7955 


0.125 


0.017265 


43/30 


1.1787 


1 


0.245065 


3 


3.9210 











TABLE I: Coefficients of the response function Eq.JSJ and the 
energy Eq.@. Results are accurate up to six digits. The coef- 
ficient p is a rational function of 7. In * (7 = 0.375)a rational 
number for p accurate to six digits could not be guessed. 



calculation is given in Appendix lAl The final result is 
1 

_ 2 2 -^r(i- 7 /2) PaPp 
A c r( 7 /2) \p\ 2 -f + 



Ac 

+ om% (9) 

The coefficients 77(7) and p{j) depend only on the po- 
tential. In Table [I] some values of the coefficients for a 
range of potentials with < 7 < 2 are listed. 



B. Continuum free energy 



e-e = e -y U K 

2 t£,\\R(l) + u(l)- R(V)-u(l')P 
1 \ 



\R(l)-R(V) 



(7) 



where u(l) is a small displacement of the particle 1 in the 
plane of the surface from its equilibrium position i?(l), 
and therefore a tangent vector to the surface JC. The 
elastic tensor n a ( g(l, 1') is defined as the leading term in 
an expansion of Eq. J2J , 



E-E = ^-J2 n ^(l, I')«a(l)«/J(l') 



(8) 



i,i' 



In deriving Eq.Q), we assume a constraint of fixed area 
per particle, enforced by a uniform background charge 
density or boundary conditions. This eliminates the term 
linear in u a (l). The physical properties of response func- 
tions are better studied in Fourier space. The detailed 



When the deviations from the ground state are small, 
the long wavelength lattice deformations may be de- 
scribed by a continuous Landau elastic energy 



F(u) = J 



d 2 r 



fiu 2 a/3 



A 2 
2 Uaa 



(10) 



The couplings A and /i are the usual Lame coefficients. 
The strain tensor u a p is defined by 



(11) 



where w(x) are the small displacements of Eq.JSJ). The 
elastic tensor Eq.©, within Landau elastic theory, is 
then 



e 2 U a /3(p) = Ac(fJ,\p\ 2 5 a /3 + (A + n)p a Pf3) 



(12) 



A comparison with Eq. immediately yields an explicit 
expression for the elastic constants of the crystal 
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1+7/2 ' 
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(14) 
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FIG. 1: (Color online:) Construction of an (n,m) icosadelta- 
hedral lattice. The filled circles indicate two nearest-neighbor 
five fold disclinations. Because these defects sit on the vertices 
of an icosahedron, they are separated by a geodesic distance 
7?cos _1 (l/\/5), where R is the sphere radius. 
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where Y is Young's modulus. The result A = oo is equiva- 
lent to a divergent compressional sound velocity as p — > 
and for 7 = 1 is just a statement of the incompressibil- 
ity of a two-dimensional Wigner crystal. Alternatively 
we can allow for wavevector-dependent elastic constants 
fj,(p) and X(p) in Ea. ((T2")l . In this case X(p) diverges as 

0, \(p) - 22 ' T - r(1 ^ /2) 1 



P 



given by Eq. lfKfy . 



A C r( 7 /2) pz 



while lim p _>oM(p) is 



FIG. 2: (Color online:) Difference in energy of (n, n) and (n, 0) 
configurations. As the number of particles in the two config- 
urations is always different (at least for relatively small n), 
we fitted the energy dependence on the number of particles 
for the two configurations, and then we computed the energy 
difference from the fitting curves. The energies are plotted 
in the inset to give an idea of the relative scale of the energy 

difference. Results are for a power law potential with 7 = 1.5 

2 

and energies are plotted in units of -gy. 



C. Spherical Crystals 

Spherical crystals have many properties not shared by 
planar ones, one of the most remarkable being that there 
is an excess of twelve positive (five- fold) disclinations. 
These disclinations repel, and the simplest spherical crys- 
tals will be those having the minimum number of defects 
(12) located at the vertices of an icosahedron. Triangular 
lattices on the sphere with an icosahedral defect pattern 
are classified by a pair of integers (n,m), as illustrated 
in Fig. ^ The path from one disclination to a neigh- 
boring disclination for an (n, m) icosadeltahedral lattice 
consists of n straight steps, a subsequent 60° turn, and 
m final straight steps. The geodesic distance between 
nearest-neighbor disclinations on a sphere of radius R is 
d = i?cos _1 (l/\/5). The total number of particles M on 
the sphere described by this (n, m)-lattice is given by [H| 

M = 10(n 2 + m 2 + nm) + 2 . (15) 

Such (n, m) configurations are believed to be ground 
states for relatively small numbers [M < 300, say) 
of particles interacting through a Coulomb potential 
US US EB OS 03 • The energy of discrete particle ar- 
rays described by Eq.JSJ) can be evaluated by starting 
with some configuration close to an (n, m) one and re- 
laxing it to find a minimum. It is found that the (n, m) 
configurations are always local minima. Whether these 
icosahedral configurations are global minima as well will 



be analyzed later. Results for the energy E(M) are shown 
in the inset to Fig. [3 

From Fig. [21 it is clear that energies grow very fast 
for increasing volume. More interestingly, the (n, 0) and 
(n, n) configurations show a growing difference in energy 
for increasing volume size, implying that the energy of 
icosahedral configurations does not tend to a universal 
value for large numbers of particles but rather remains 
sensitive to the (n,ra) configuration, a result also noted 
by other authors |l9l l39| . Further insight comes from 
investigating the distribution of energy. Plots of the local 
electrostatic energy, the electrostatic energy at point x on 
the sphere, as defined in Eq.© are shown in Fig. 40J . 

From Fig. [3] it should be noted that the triangles ob- 
tained by the Voronoi-Delaunay construction, after min- 
imization of the potential Eq.JSJ), are very close to equi- 
lateral. 

The distribution of the local energies for the (n, 0) and 
(n, n) configurations are very different. The (n, 0) config- 
uration shows maximum energies along the paths joining 
the defects. The (n, n) configuration, on the other hand, 
has its maximum energies along the directions defined by 
the triangles formed by three nearest neighbor disclina- 
tions. The size of these regions of differing electrostatic 
energy turns out to scale with system size, making it very 
plausible that there might be small differences in the en- 
ergy per particle for (n, 0) and (n, n) configurations in 
the limit n — > 00. This point will be discussed in more 
detail in coming sections. 
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FIG. 3: (Color online:) Potential energy distribution for a 
(n, 0) configuration with n=l and M=f002 (top) and a (n, n) 
configuration with n=6 and M=f082 (bottom). 



III. THE GEOMETRIC APPROACH 



The minimization of functional forms like Eq.(|5j) is 
hampered by the computational complexity of the prob- 
lem, which is exponential in the particle number for 
spherical crystals |4l|. This difficulty, which is made 
worse by the "geometric frustration" associated with 
packing particles on the sphere, limits direct approaches 
to minimizing the energy to systems having a small num- 
ber of particles, even if much larger computer resources 
become available. 

One way to overcome these difficulties is to substan- 
tially reduce the number of degrees of freedom that need 
to be considered. An approach focusing on the topo- 
logical defects as degrees of freedom, rather than on the 
actual particles, was proposed in Refs.|2ll l4q. Some as- 
pects of this formalism are now described. The concep- 
tual issues and developments presented in this section are 
applicable to crystals in any topography. Some of the re- 
sults given here have already appeared in brief form in 
Ref.Ba. 



A. Effective free energy 

The elastic energy of a curved crystal may be obtained 
by writing in a parametrization-invariant way the results 
for a flat crystal. If the metric of the curved surface is 
g a fj (with determinant g), the energy reads 

n/T = E + j(p\±p) + ^{p\±p)+E s a 2 (s\s) , (16) 

where Eq is the energy of a defect-free monolayer, 
(A\B) = f (Px^/gAB, p(x) = K(x) - s(x) with K(x) 
the Gaussian curvature and s(x) the disclination density 
s(x) = g^= Yli=i qiS(x~Xi). Here Y is a Young's modu- 
lus and Ka is a hexatic stiffness constant. We have added 
a core energy term to account for the short-distance 
physics of disclination defects. The quantity -jjp(x) has 
the meaning 

i^(*) = /d 2 xV?(^) x , x ,p(x'). (17) 

A similar expression can be defined for ^p{x) = i/)(x). 
Here, G(x, x') = (^s) , is a shorthand notation for a 
Green's function which obeys V 2 G(x, x') = y/g~8(x — x'), 
where V 2 is the covariant Laplacian. Positive/negative 
disclinations are attracted to positive/negative curvature 
regions respectively. We note that at finite temperature, 
an additional term proportional to 

(K\^K) , (18) 

arises from the short distance behavior of the measure 
(the Liouville anomaly) ^(| . This term can be safely ig- 
nored in the present analysis which focuses on zero tem- 
perature. 

The defect part of the free energy Eq. (|16|) will be used 
in a simplified form in the crystalline phase. In that 
phase the hexatic term can be incorporated into a core 
energy contribution proportional to the total number of 
defects. The energy we need to minimize becomes 

E = E + ^(p\^p) + NE c , (19) 

where N is the total number of disclinations of core en- 
ergy E c . 

If the disclination density were continuous, instead of 
being composed of discrete objects, configurations of de- 
fects such that 

p = =S> s(x) = K (x) , (20) 

would be absolute minima of the free energy Ea. (|19fl . 
In general, defects tend to arrange themselves on curved 
surfaces to screen the Gaussian curvature as efficiently 
as possible consistent with their discrete topological 
charges. 

The free energy just discussed can also be applied to 
fluctuating geometries, as in the case of fluid or hexatic 
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membranes (see El for reviews) . If Young's 

modulus vanishes, corresponding to a proliferation of un- 
bound dislocations, one obtains the free energy of an hex- 
atic membrane Ho, IH3 . 



IV. 



GEOMETRIC FORMALISM ON THE 
SPHERE 



Spherical substrates provide the simplest example of 
the problem of crystals on curved surfaces. The study of 
spherical crystals is simplified by two important proper- 
ties: there is a unique scale with dimensions of length, 
the radius R, and there is a fixed excess disclinicity of 
twelve following from the Gauss-Bonnet theorem 



N 

/ d 2 x >/fl(x)s(x) = 4tt -»■ ft = 12 



(21) 



The free energy Eq.JTSJ, applied to the sphere, 
is tractable analytically because the inverse square- 
Laplacian operator on a sphere of radius R can be com- 
puted explicitly. It is shown in (2l| that the Green's func- 
tion for the square Laplacian, in spherical coordinates 
(0, </>), has the following simple form on a unit sphere: 



x(0°,^W) = i + 



dz 



In ; 



1 



(22) 



where (3 is the geodesic distance between two disclina- 
tions located at (6 a ,(j) a ) and (6» b > & ), 

cos [3 = cos 6 a cos 6> h + sin a sin 6» b cos(0 a - <j> b ) . (23) 

The total energy of a spherical crystal with an arbi- 
trary number of disclinations follows from Ea. (|19(l and 
Ea.p5)l and has the simple form [2l[ 



ttY 



N N 



2E(Y) =E Q + —R 2 ^ ^ WXW^V ,^) + NE c , 

i—1 j — 1 

(24) 

where 0i}j=i,... ,j\r are the coordinates of AT defects 
and we restrict ourselves to 5-fold (g, = +l)and 7- fold 
(<ft = —1) defects. The quantity Eq is the zero point 
energy and is defined in Eq. |JB}. Although 5 and 7- 
fold disclinations will in general have different core en- 
ergies }56j . we assume equal core energies here for sim- 
plicity. What matters for our calculations in any case 
is the dislocation core energy E c i, which we take to be 
Ed = E$ + £7 = 2E C . 

The value of the Young's modulus and the flat space 
ground state energy £0 have been computed in Sect. Ill Al 
When the sphere radius R is large compared to the par- 
ticle spacing a, we can use flat space values of Y and the 
flat space energy Eq(M) associated with a finite num- 
ber of particles M. To obtain the leading terms in the 
expansion of the ground state energy for large but finite 
M, the precise compactification of the plane employed 



is irrelevant - it may be achieved by periodic boundary 
conditions, for example. For a sufficiently large plane 
the finite size effects will be negligible. The density a 
of particles is then M divided by the total surface of the 
compact plane, taken to be the surface area of the sphere 
of radius i?, 



1/Ac 



AI 

~s 



S = AnR 2 



(25) 



From Ea. i|13|) the expression for the Young's modulus 
suitable for M particles on a spherical crystal of radius 
R with < 7 < 2 is then 



Y = 4/x 



477(7)M 1+ ^/ 2 e 2 

( 47r )l+ 7 /2 R 2+j 



(26) 



One remaining detail is the divergent contribution to the 
energy Eq in Eq.©. Since the divergent part comes 
solely from the zero mode, the spatial variations in the 
density of the actual distribution are irrelevant. It may 
therefore be computed for a uniform density of charges. 
The divergent part is identical to the energy of a constant 
continuum of charges as described by the density Ea. (|25|l . 
We now evaluate this divergent part of the energy on a 
sphere, instead of a plane. 



E D 



Me' 



7 /2 



2 2- 7 



V^x)p(x)^— _p( x ')Vff(x') (27) 
M 2 



|X - x'lT 

o2 



27-1(2 - 7) Ri 



The divergent part has thus been regularized, and the 
energy is finite and well-defined for all M < 00. 

Note that for the case 7 < 2 of primary interest to us 
here, E D - M 2 ~i 'l 2 (M / 'S) 1 ,/2 . Hence E D is not simply a 
function of the particle density M / S, as one would expect 
for a short range interaction. 



A. The Energy of Spherical Crystals 

Upon substituting the elastic constant of Eq. (|26|l into 
Ea. (|24|) . one arrives at 



ttY 



N N 



2E = E + —R^^qiq^,^;^ + NE C 



i=i j=i 



(28) 

where Eq is defined in Eq. and the function 

C(i\ - ■ ■iff) depends on the position = (#i,0i) etc. 
of the A^ disclination charges and is universal with re- 
spect to the potential. The total energy of a spherical 
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crystal, including the contributions to Eq is then 



2E> 



TOT 



(M) = 



M 2 



27-1(2 - 7 ) 

477(7) 7r 
(4tt) 1 +^/ 2 36 

NE r . 



0{l) 
(4tt)7/2 



C(h, ■ ■■ ,in) 



M l+7/2 




Note that the leading correction to the zero mode energy 
proportional to M varies as M 1+7 / 2 , and depends both 
on the flat space function #(7) and on the C-coefncient 



N N 



C(ii 



,*aO = EE ^x(#\^;^ j ) 



(30) 



3=1 i=i 



associated with a particular configuration of disclina- 
tions. 

Note that the core energies contribute to the second 
sub-leading coefficient. For short-range potentials, such 
as 7 > 2, the ground energy is extensive, and the leading 
term varies as M 1+7 / 2 . 

The extensive nature of the M 1+7 / 2 term becomes 
clear upon noting that 



FIG. 4: (Color Online:) Illustration of the calculation done in 
the text. The energy of the discrete (n, n) configuration on the 
left is extrapolated for large M and compared to the energy 
computed with the continuum model on the right. While in 
the continuum model only twelve degrees of freedom (the 12 
disclinations) need to be considered, the direct calculation of 
a family of discrete models requires the consideration of the 
full lattice and a careful extrapolation of the energies to large 
M. 



a frozen topography we expect an expansion along the 
lines of Ea.pty. 



2B 



tot(M) = a M 2 - atW** 1 -* 



(34) 



Ri a 7 + 2 



(31) 



where a is the particle spacing. Comparison with Eq. (2H 
shows that the dimension of Young's modulus Y arises 
solely from the lattice constant a and the electric charge 
e, consistent with elastic constants arising from physics 
on the scale of the lattice constant in an essentially flat 
geometry. This observation is now generalized to the rest 
of the couplings discussed in the previous section. 
For the hexatic term, in Ea. (|16|l . we have 



K A 1 
— (P\^P) 



K A R° 



e 2 1 p 2 



CP 



(32) 



Since core energies arise as short-distance divergence's 
similar to the hexatic term, they are a sub-leading con- 
tribution. For a fluid membrane not on a frozen topogra- 
phy, Helfrich terms arising from the extrinsic curvature 
H (x) as well as the Gaussian curvature can be important. 
These scale in a way similar to the hexatic term, 



k j dxy/gH 2 
k g J dxyfgK(x) 



K [ dxJTjX; = kR° - Mil 2 ^- , 
J vy i? 2 Rf 

K f cWg-^ = kR° - Mi' 2 — 
J v y R 2 Ri 



(33) 



Both terms would therefore contribute to the same order 
in the M expansion as the hexatic term, although the 
last term is purely topological. For crystals embedded in 



The nonextensive term a$M 2 arises from the long range 
interactions. The next extensive contribution comes from 
the interaction between Gaussian curvature and defects 
as well as the extensive energy per particle in flat space. 
Hexatic terms and bending rigidity contributions are 
higher order in 1/M and can be absorbed into a re- 
definition of the disclination core energy. Core energies 
also depend on non- universal details of the short-distance 
physics. Core energies are included explicitly in Eas. i|19|) 
and I2IJ) ■ 

The results presented so far are strictly for systems at 
zero temperature. In systems with short range interac- 
tions, the elastic constants can be strongly temperature- 
dependent. An extreme example is hard disks of radius 
ao, which may be viewed as a limiting case of a power law 
potential of the form V{r) ~ eo (ao/r) 7 , with 7 — > 00. In 
this case, the elastic constants are strictly proportional 
to temperature. It is straightforward, however, to adapt 
the techniques of this paper to the simpler problem of 
short range pair potentials. 



B. Energies of Icosahedral configurations 

The configuration on a sphere with the minimum num- 
ber of charge ±1 defects is twelve +1 (5-valent) disclina- 
tions, which minimize their energy by sitting at the ver- 
tices of an icosahedron y. The energies of such config- 
urations will be computed for the discrete spherical tes- 
sellations described in Sect. Ill Gl and compared with the 
predictions of continuum elastic theory, as illustrated in 
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Fig-El It is well established that for sufficiently large val- 
ues of M configurations with more than 12 disclinations 
(i.e., those with "grain boundary scars") have lower ener- 
gies 0,0, El- It is of interest, however, to study sim- 
ple icosahedral configurations for large M, as metastablc 
states with a well defined energy. 

Within the continuum elastic theory it can be shown 
that twelve disclinations at the vertices of an icosahedron 
minimize the energy plj when no further defects are al- 
lowed. The C-coefficient of Eq . (|29l) for this configuration 
of defects has been computed in |21| [64j 



C(y) = 0.6043 . 



(35) 



y here stands for a particle configuration with 12 defects 
at the vertices of an icosahedron. Using the energy of 
Ea. (|29fl . the coefficient 01(7, y) appearing in the expan- 
sion of Eq. (|34|l may be computed, with the results shown 
in Fig. EJand Table EU 

From the results described in Sect. Ill 01 the a\ coeffi- 
cient may be extrapolated to very large numbers of par- 
ticles using the expansion derived from Ea. (|34|l . Indeed, 
as shown in Fig. El plots of 



i{M) 



_ 2R'<E TOT {M)/e 2 - a (7)Af 2 



(36) 



vs 1/M are linear, with a slope that determines 01(7) 
and an intercept related to the higher order core energy- 
like contribution. The results of these extrapolations are 
shown in Table [H] The agreement between the contin- 
uum elastic theory and the explicit computation for the 
(n, n) configuration is remarkable, holding to almost five 
significant figures. For the (n, 0) lattice there is agree- 
ment to four significant figures. This agreement is even 
more striking when it is recalled that the ai-coefficient 
is obtained after subtraction of the term oo(7)M 2 , as 
illustrated in Fig. |21 Furthermore, in the range from 
7 = 0.125 to 7 = 1.875, all the significant digits vary and 
yet the accuracy of the calculation is virtually indepen- 
dent of 7. 



4.9 



3.9 



?f 2.9 



1.9 



0.9 



0.5 1 1.5 2 

J 



3.0e-04 



2.0e-04 



1.0e-04 



0.0e+00 




FIG. 5: (Color online:) Energy coefficient a\ as a func- 
tion of gamma (solid line) and from the numerical results 
with (n, m) configurations (filled circles), for the icosahe- 
dral configurations. Plot of 01(7, y) - ai(j,y) l ~ n ' n ' > (cir- 
cles), 01(7,3;) - ai(7,J) (n ' 0) (diamonds), ai(-y,y) {n ' n) - 
ai(7,3 7 ) (n '° ) (squares). 



V. THOMSON PROBLEM WITH A 
CONTINUOUS DISTRIBUTION OF 
DISLOCATIONS 



C. The Energy difference of the (n, m) lattices 

The ai-coefficient computed within our continuum 
elastic approach above does not depend on the icosadelta- 
hedral class (n, m). Results from the direct minimization 
of particles do, however, show a weak dependence (in the 
Ath significant digit) on the particular (n, m) configura- 
tion, as is apparent from Fig. |5] and Table [D] It should 
be noted that the discrepancy from the continuum re- 
sult has a well defined sign, and is therefore reasonably 
attributed to a term not present in the energy expansion. 



When the number of particles is extremely large, the 
minimum energy configurations can be approximated by 
a closed analytical form, upon assuming a continuous dis- 
tribution of defects. Only the sphere will be worked out 
here, but other curved surfaces can be treated in a very 
similar fashion. 

The formal elimination of the geometric frustration in- 
troduced by the Gaussian curvature may be formulated 
as a concrete set of equations in the case of the sphere. 
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1/M 



0.004 



FIG. 6: (Color online:) Numerical estimate of e(M) as a 
function of 1/M for (n, 0) and (n, n) icosadeltahedral lattices 
with 7 = (1.5,0.5). 
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(n,n) 


(n,0) 


1.875 


4.45118 


4.45110(4) 


4.45095(4) 


1.75 


2.47175 


2.47166(3) 


2.47150(3) 


1.625 


1.82629 


1.82621(2) 


1.82603(2) 


1.5 


1.51473 


1.51454(2) 


1.51445(2) 


1.375 


1.33695 


1.33683(4) 


1.33667(4) 


1.25 


1.22617 


1.22599(7) 


1.22589(7) 


1.125 


1.15366 


1.1535(2) 


1.15340(2) 


1.0 


1.10494 


1.10482(3) 


1.10464(3) 


0.875 


1.07187 


1.07174(3) 


1.07160(3) 


0.75 


1.04940 


1.04921(6) 


1.04910(6) 


0.625 


1.03421 


1.03413(5) 


1.03398(5) 


0.5 


1.02392 


1.02390(4) 


1.02372(4) 


0.375 


1.01672 


1.01663(6) 


1.01656(6) 


0.25 


1.01115 


1.01106(3) 


1.01103(3) 


0.125 


1.00595 


1.00592(2) 


1.00589(2) 



TABLE II: Numerical values of the coefficient 01(7, y) (twelve 
disclinations on the vertices of an icosahedron) using the C- 
coefficient from Ea. l|35^ . The same coefficients from the (n, n) 
and (n, 0) lattices. 



We shall use the identity 

N 



(37) 



n 00 l N 

1=1 m=—l i=l 
00 I N 

^« + ^E E ^T(M)E^(^): 



1=1 m——l i—1 

which follows from the topological constraint Ea. (l21ll . 



7 ai(7>0) 7 ai(l,G) 7 ai(l,S) 7 ai(l,G) 



1.875 


4.45227 


0.875 


1.07297 


1.75 


2.47289 


0.75 


1.05044 


1.625 


1.82746 


0.625 


1.03515 


1.5 


1.51592 


0.5 


1.02473 


1.375 


1.33815 


0.375 


1.01737 


1.25 


1.22737 


0.25 


1.01161 


1.125 


1.15485 


0.125 


1.00620 


1 


1.10610 







TABLE III: Value of the ai coefficients for the Q configuration 
Eq.lEHJ. 



Provided a disclination configuration exists such that 

N 



5>il(^) = o 



(38) 



for each {I > l,m), the disclination density completely 
screens the Gaussian curvature. A configuration of de- 
fects satisfying Ea. (|38|l is an absolute minimum of the 
elastic energy, a result easily understood by writing the 
energy in the form 



E — En 



ir 2 Y 
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-R 



EE 

l—l m——l 



\T N a-Y l 



■NE B 



where the zero point energy En in Ea. (|16fl is kept. A 
configuration satisfying Eq. (|38|) will be denoted by Q. 
For this hypothetical configuration, the C-coefficient in 
Ea. (|29|) vanishes, although there is now a large contri- 
bution (linear in R) from the dislocation core energies 
represented by the last term of Eq. (|39|l . 

The Q configuration may be characterized more explic- 
itly. It consists of a density of dislocations that converges 
to the local Gaussian curvature. It can be shown that 
upon approximating the dislocations (each regarded as a 
disclination dipole with spacing a) as a continuum dis- 
tribution, this dislocation density for a sphere becomes 



1 6 

b (^) = ^E cot M^)R • 



(40) 



fe=i 



The summation here runs over the six coordinates of 
the northern hemisphere of an icosahedron ((0, 0) and 
[dy, ^p), where 6y — arccos(^) and otk is the angle 9 
relative to a coordinate system with the north pole lo- 
cated at (0y, 2 ^ L ) for k = 1, • • • ,5. This angle is given 
implicitly by 

2ir 

cos[afe(6', ip)] = cos(6>) cos(9y)—sm(8) sm(9y) cos( — k+ip) 

5 

( 41 ) 

The implicit form of Eq. I|4U|I can be further simplified 



2tt 

— sin(6>v) sin(— k + ip)eg 
5 



(42) 



2vr. 



+ [{cos(6*3;) sin# + sin(0j;)} cos9cos(ip + — k)]e v 





10 



where f k (9, <p) 



-i(a k (6,ip)) 



Close to one of the 12 disclinations with charges = 
Eq.l^DJ predicts a singularity in the dislocation density 
451 



(43) 



2vri?a 



For small angles, close to each disclination, there is a 
short-distance singularity 



^ = we 



(44) 



in agreement with known results in flat space. 

Ea. H40|) represents a continuous distribution of dislo- 
cations, and neglects both dislocation discreteness and 
their mutual interactions. It represents six families of 
dislocations with azimuthal Burgers' vectors associated 
with antipodal pairs of the 12 original disclinations in 
the icosahedron. When discreteness and interactions are 
taken in account, we expect these dislocations to con- 
dense into grain boundary arms, containing with quan- 
tized Burgers' vectors and variable spacing in the radial 
direction |2j Hj, |53 • The total number of discrete arms 
remains, therefore, the variable that needs to be deter- 
mined for a discrete solution of the Thomson problem. 




- 



FIG. 7: Results of a minimization of 500 particles interacting 
with a Coulomb potential, showing the appearance of scars. 



A. The intermediate regime 

Within the continuum elastic approach, the dominant 
configurations for a small number of particles are 12 de- 
fects with an icosahedral symmetry [21j. We have just 
seen, however, that adding a continuous distribution of 
dislocations, as might be appropriate when the particle 
number is large, can more efficiently screen the Gaus- 
sian curvature on a sphere. The natural problem then 
becomes to determine the precise structure of the de- 
fect arrays for intermediate numbers of particles when 
the discreteness of interacting dislocations is taken into 
account. 

We note first that the particular arrangement of defects 
dominating in this regime will not be fully universal. The 
particular array structure favored can vary from system 
to system with fixed particle number, depending, e.g., 
on details such as the dislocation core energy. This re- 
sult may be traced back to the M-expansion of Eq. (|34ll . 
in which the sub-leading terms which depend on non- 
universal properties influence the dominant terms for fi- 
nite values of M. Some typical defect configurations ob- 
tained by direct minimization of particles on the sphere 
are shown in Fig. and show incipient scars, already 
at number of particles of 500 (in 0, the minimum 
number of particles where scars are systematically found 
is predicted around 400). By using the geometrical model 
described in this paper, where the energy is parameter- 
ized just b y a Young's modulus and a dislocation core 
energy |2ll l22j one can simulate larger particle numbers 
and one obtains results as in Fig. EI Note the occurrence 



of low energy configurations with scars (to = 2) in one in- 
stance and pentagonal buttons (to = 5) in another. The 
dislocation spacing decreases the further a dislocation is 
from the central disclination. 

An overview of previous results involving grain bound- 
ary scars is presented in Fig. If a disclination is placed 
on a perfect crystal, no additional defects will appear if 
the disclination is located on the tip of a cone with total 
Gaussian curvature equal to the disclination charge. If a 
disclination is forced into a flat monolayer, then to low 
angle grain boundaries, with constant spacing between 
dislocations as shown in Fig. [5] and grains goi ng all the 
way to the boundary, will be favored (see [44J for a de- 
tailed discussion). In the intermediate situation where 
a finite Gaussian curvature is spread over a finite area, 
as in the case of a spherical cap, a disclination arises at 
the center of the cap, and finite length grain boundaries 
stretched out over an area of ^R 2 with variable spacing 
dominate, again as illustrated in Fig. Since several 
non-universal features, related to the size of the core en- 
ergies, commensurability properties and so on, will have 
an important effect in this regime of M, the previous ar- 
gument should describe the general trends and will be 
realized in an approximate form only. 

Additional results may be obtained for the number of 
arms within the grain boundary, the actual variable spac- 
ing between dislocations within the grain and the length 
of the grains as a function of the number of particles. 
The detailed study of these questions will be reported 
elsewhere. 
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FIG. 8: Ground state configurations for M « 2000 particle 
obtained from the continuum elastic formalism. In the to] 
figure one finds scars (m — 2) and in the bottom pentagona 
buttons (m = 5) forming a rhombic tricontahedron. 



Crystal with no defects 




Spherical cap 



m grain boundaries 

Constant spacing of dislocations within a 




m grain boundaries 

Grain boundanes have finite length 
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When grain boundary scars appear, we can estimab 
the number of excess dislocations which decorate each of 
the 12 curvature-induced disclinations on the sphere us- 
ing ideas from Ref.[2]]]. This estimate is in reasonable 
agreement with experiments probing equilibrated assem- 
blies of polystyrene beads on water droplets j^. Con- 
sider the region surrounding one of the 12 excess discli- 
nations, with charge s — 27r/6, centered on the north 
pole. As discussed in Ref.[21|. we expect the stresses 
and strains at a fixed geodesic distance r from the pole 
on a sphere of radius R to be controlled by an effective 
disclination charge 



FIG. 9: Schematic illustrating the genesis of grain boundary 
scars. A disclination is first constructed from a perfect lattice. 
If this disclination is placed on a tip of a cone, with a delta 
function of Gaussian curvature balancing the defect charge, 
then no additional defects form. If the crystal is forced into 
a monolayer, grain boundaries radiating out of the disclina- 
tion radiate all the way to the boundary. In the intermedi- 
ate regime of constant non-zero Gaussian curvature, m grain 
boundaries of finite length and variable spacing of dislocations 
form. 



s ef f(r) = s- d<W &r'^(r)K 



7 r/3-4 7 rsin 2 (— ). 



(45) 



Here the Gaussian curvature is K = 1/i? 2 and the metric 
tensor associated with spherical polar coordinates (r, <p) , 
with distance element ds 2 = d 2 r + R 2 sin 2 (r/R)d 2 cj), gives 
\J g(r) = R sin (r/R). Suppose m grain boundaries radi- 
ate from the disclination at the north pole. Then, in an 
approximation which neglects interactions between the 
individual arms, the spacing between the dislocations in 
these grains is [2l| 



which implies an effective dislocation density 



n d {r) 



1 
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ma 
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2tt 






cos 


ma 





47rsin 2 (r/2i?) 



r 
R 



This density vanishes when r — > r c , where 
r c = Rcos- 1 5/6 » i?(33.56°) , 



(47) 



(48) 



(46) 



which is the distance at which the m grain boundaries ter- 
minate. The total number of dislocations residing within 
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this radius is thus 



VI. CONCLUSIONS 
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sin 2 (r/2i?)dr 



/ TT- 5 cos" 1 (5/6) 
0A08(R/a) . 



(R/a) 



(49) 



As discussed in Ref . j37l| . it is also of interest to consider 
2tt disclination defects (appropriate to crystals of tilted 
molecules 59]) on the sphere. The icosahedral configura- 
tion of 12 s — 2tt/6 disclinations is now replaced by just 
two s = 2-7T disclinations at the north and south poles. 
Using the approximation discussed above, it is straight- 
forward to show that the density of dislocations in each 
of m (noninteracting) grain boundary arms now reads 



n d (r) 



2?r ,r , 

cos ( — ) 

ma R 



(50) 



This density vanishes at r c = corresponding to a 
hemisphere of area on the sphere for each cluster of arms. 

It is of considerable interest to repeat the above cal- 
culation for a square lattice, as found for example in the 
protein surface layers (s-layers) of some bacteria |47ll48j. 
In this case the basic disclination has s — 27r/4. The 
effective dislocation density becomes 



n d {r) 



1 

W) 

1 

ma 
2tt 
ma 



- -47rsin 2 (r/2i?) 



cos 



R 



3/4 



(51) 



This density vanishes when r — > r c , where 

r c = Rcos- 1 3/4 » i?(41.4°) , (52) 

which is the distance at which the m grain boundaries 
terminate. The longer angular length of square lattice 
scars reflects the larger initial disclination charge (90°) 
that must be screened. The total number of dislocations 
residing within this radius is thus 
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n d {r)dr 



47T 



sin 2 (r/2i?)dr 



3 cos" 1 (3/4) 



(R/a) 



0.75(R/c 



(53) 



Thus the angular length of scars and the total number of 
excess dislocations is a measure of the underlying topol- 
ogy of the lattice tiling the sphere. 



A. Summary of Results 

In this section we summarize the most relevant results 
obtained from the analysis presented earlier. 

In Sect. ^ we treated several properties of planar and 
spherical crystals which were subsequently used to test 
our continuum elastic formalism. We computed the en- 
ergy Eq.© and elastic tensor Eq.® for triangular lat- 
tices in flat space for a general long range power law po- 
tential of the type Eq.J5J. The continuum elastic formal- 
ism, where defects such as disclinations and disclination 
dipoles = dislocations are the relevant degrees of freedom 
and six-coordinated particles are treated as a continuous 
elastic background, was discussed in Sect. IIHI It was 
shown that the total energy is expressible as an expan- 
sion in powers in the total number of particles [Ea. 1)29(1 ]: 



2E 



M 2 



2^- 1 (2- 7 ) 
a 2 ( 1 \{q i }^ 1 ,, N )IVp/ 2 



a 1 ( 7 |{g I } l =i,, J v)M 1+ ^ 2 (54) 



e 



where each coefficient has a clear geometric interpreta- 
tion in terms of continuum results. 

Our approach was illustrated for the generalized 
Thomson problem in Sect IIVI Using the elastic con- 
stants computed in flat space, the continuum elastic the- 
ory gives concrete energy predictions, with no fitting pa- 
rameters, as a series expansion in the total number of 
particles M, which can be compared with the energies ob- 
tained numerically for spherical crystals. We find agree- 
ment to 5 significant figures for (n, n) icosadeltahedral 
lattices and to 4 significant figures for (n, 0) lattices, as 
presented in Table [H] Only a small discrepancy, of or- 
der the difference between (n, 0) and (n, n) tessellations, 
separates the continuum results from results for actual 
spherical crystals. 

The limit of a very large number of particles M was 
dealt with in Sect[V] A "Debye-Huckel" type formulation 
where dislocations are treated in a smeared continuum 
density of Burgers' vectors was proposed. In Ref. (5?| an 
explicit solution for the actual defect distribution without 
assuming a continuum of dislocations was proposed and 
it was shown that certain dislocation grain boundaries 
have a C-coefficient that vanishes in the limit of a very 
large number of particles (R/a — > oo). This solution is 
a discrete version of the continuum solution presented in 
this paper, and incorporates the discreteness of the dis- 
location positions and charges and their mutual interac- 
tions. We should mention that an alternative scenario for 
the Thomson problem has been proposed 58], where at 
some finite value of number of particles an instability to 
a "spontaneously magnetized" state is predicted. Based 
on the results presented in this paper and in Ref . [3(| |57j 
we conclude that such instability does not appear for 
the generalized Thomson problem. It is possible that an 
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instability of the type predicted in [58| may appear for 
charges on spheres under other type of constraints. 

The intermediate regime was discussed in the last sec- 
tion and it was shown that the underlying universality of 
the result competes with several non-universal features 
of the problem. 



B. Outlook 

The main goal of this paper was to introduce a con- 
tinuum elastic approach to address the problem of two- 
dimensional crystals in frozen topographies. The formal- 
ism has been explicitly applied to the sphere, but it ap- 
pears general enough to be applicable to a variety of other 
geometries. The case of crystalline order on a torus is 
currently under exploration. 

We hope this presentation will inspire further work on 
the problem of crystals on curved topographies. The long 
range pair interactions on a sphere studied here certainly 
do not exhaust the possible problems. 
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APPENDIX A: THE EVALUATION OF 
POTENTIALS 

The details of the computation of the energy and the 
elastic tensor, Eqs.(0 and Eq.© are described in detail. 
The approach followed is a gene ralization of the one used 
by Bonsall and Maradudin |60| . See also ref [6l| . 



1. Computation of the energy 

The energy for a system of M particles located at po- 
sitions R(V) of a 2d Bravais lattice, defined by vectors an 
and in direct and reciprocal space 



R(l) = hai + l 2 &2 
G(h) = hiei + h 2 e 2 



(Al) 
(A2) 



is given by Eq.JSJ 



2 M 1 

E = 6 - V 

M%- lira Y \ = ^-EM . (A3) 

2 3~ofe\g- m t 2 



To efficiently perform the sum (a generalization of the 
Ewald method), we separate short and long distances 
contributions, since they give rise to different singular 
behavior. This may be achieved by the identity Ea. ljBlfl 



'2\ 
2> <t 



\x-R{i)\i r( 

+ f 'dtt-i+7/Sg-tiif-fia 
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CT 7/2 



r p 7/a _i((7|a-j2(i)| a ) 



r 7/2 r 1 



r(i) J 



(A4) 



The definition of the Misra functions (p n is given below 
in Eq. (Bj^- 
Using the Poisson summation formula Eq. i|B5|l , the last 
term in Eq. I|A4(I can also be expressed in terms of Misra 
functions, 



fl(i)l 2 



(A5) 



Ac<? Jl 
G 



Ana- 



e P-7/a(-^r) 



G 



where G are the vectors in reciprocal space of the Bravais 
lattice, and Ac is the area of the unit cell. 

Upon combining Ea. l|A4(l and Eq. l|A5|l . the energy 
Ea. (|A3p becomes 



r(i) 



— J2^/2-i(a\x-R(l)\ 2 ) 

1^=0 

\G\' 



7 r( 7 /2) 



A c r( 7 /2) 



e V-7/a(-^r) 



(A6) 



Although the limit x — > is in general convergent, the 
term G — requires special attention. This term has to 
be treated separately by considering \G\ small but non- 
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vanishing, 



we a 



2 ,t7/2-1 



iGx 



r(7/2) 



V-7/a(-^r) 



7r2 2 -^r(i- 7 /2) c2 
|G|2-7r( 7 /2) 

A c r( 7 /2)(i- 7 /2) e 

+ 0(\G\). (A7) 

Thus, the singularity as G — > 0, associated with the large 
distance behavior, has been explicitly isolated. 

The results derived so far are completely general, valid 
for any Bravais lattice. Since only the triangular lattice 
is relevant to this paper, complete results for other Bra- 
vais lattices will be published elsewhere. The two vectors 
ai , a 2 defining the triangular lattice in direct space, and 
the two vectors ei,e2 in reciprocal space are taken as 



ai = (a,0) , ex = — (1, — — ) 
a 6 



a 2 



„a \/3a . 



2vr , 

e 2 = — 0,2^- 
a 3 



(A8) 



where a is the lattice constant. Further simplification is 
achieved by choosing a as a — -j^ (where Ac — -^a 2 ), 
so that the argument in the Misra functions has a similar 
form for the sums in both direct and reciprocal space. 
Thus, 



The property R(l) - R(Y) = R(l - Y) - R(0) implies 
translational invariance 11(1,1') = 11(1 — l',0). The re- 
sponse function is better studied in Fourier space. The 
Fourier transformed elastic tensor can be computed from 
the identity 

n Q/3 (p) = -(S a p(d - s a/3 {6)) , (Ai3) 

with S a p defined as 



S a p(p) 



x^o d a dp 
lim 



x-R(l)\' 



x^o d a dj3 



F(%,P) = E 



s ip-(x-R(l)) 

x-R(l)\-> 



(A14) 



The function T can be computed by further using 
Eci. (|A4(l . Ea.l|A5|l and Ea. (|B5|) . with essentially the same 
steps as in previous computations, leading to the expres- 
sion 



ai? 2 (l) = -^(/fa? 
G\h) 



2/iZ 2 aia 2 + 1%b%) 



, (h\e{ + 2hih 2 e 1 e 2 - 
4cr 4ct 

= -r-{h\s? l -2hih 2 a.\a.2 

Our final form for the energy is then 
,2 r „ 4 



h\el) 



(A9) 
(A10) 



E(i) 



r(7/2) 



(— V' 2 \- 

{ A C > S(2-7) 



JL^(i))-E ^/2-i(^(i))}l 

Ac U c J 



vr 2 2 -7r(i- 7 /2) 

|gho A c |G| 2 -7T( 7 /2)' 



lim 



(AH) 



The summation over the Misra functions is exponentially 
convergent; just a few terms give a very accurate result. 



r 7/2 



r(7/2) 



E< 



-ipR(\) 



ip 7/ 2-i^\x-R(l)\ 2 ) 



^-7/2 /■! 



r(7/2) 



n 



A c F( 7 /2) E 



4a 



7/2-1 



A c r( 7 /2)^ 2( ^^ } ' 



(A15) 



Upon inserting the derivatives of the Misra functions ob- 
tained from Eq. [|B3f> . Eq.|jB4} an( i using Eq. (|A14|l . we 

have 



cl TTCrT/ 2 " 1 -., _ -> ,|p+G| 2 , 
SaMW = ~ ^ cr ( 7 / 2 ) 2^^ +G ^^ + G ^ y~7/2( 4^ ) 



2. Computation of the elastic tensor 

The response function is defined in Eq. iJSJ . Using some 
simple algebra, the explicit form for this tensor response 
function II follows 

( 2 x (fi(i)-fl(i , )) a _(fl(i)--R(r))3 
rw 7 |(fl(i)-fl(i')|^+ 4 

n^(i,i') = <! + 7 ^ 



Kfl(i)-ij(i')i^+ 2 
Ei^ n Q/3 (i,i') 



i^i' 
i = r 

(A12) 



2(7 7/2+l 
4^2+2 

r(772T 



E- 



-ip-R(\) 



^,2^\R{\)\ 2 ) 



J2R(l) a R(Y)p ^ 1+7/2 (a|i?(l)| 2 ) 



4 g^ 2+1 
2 + aT( 7 /2)' 



The full expression for the response function is then 



(A16) 
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n a/ 3(p) 



■ Z> + g )^p + G W-7/2( Tn ) ~ E G « G /3 



A c r( 7 /2) 



2(T T/2+l 



-ip-R(X) 



Aa ' A c r( 7 /2 

4c7 7/2+2 



G 

ip-fl(l) 



1) ^ 7/2 (<7|ii(l)| 2 ) - ^-^-^(e-^W - 1)^(1)^(1)^^(^1^(1)1^17) 



r 



As in the computation of the energy, it is convenient to to non-analyticities. A Taylor expansion for the G — ► 
isolate the G — > contribution since it usually gives raise contribution leads to a final expression 



22-^(1-7/2)^ 22-^(1-7/2)^ / 1 



A c r( 7 /2) bi 2 -7 



7/2.- /-P \ _ 1 

A c r( 7 /2) |pl 2 -7 Vr(l- 7 /2) { 4a> ^- l/2[ 4a' 



KG 



7/2-1 



^ c r( 7 /2) 



\P + G\' 



) 



A c T( 7 /2) 



(A18) 



Oo-7/2+i ^ 40-7/2+2 

■J2[l-oo S (p-R(l))]^ /2 (a\R(l)\ 2 ) ' ' 



r(7/2) 



MO 



r(7/2) 



£[1 - cos(p • R(l))}R(l) a R(l) p <p 1+ ^ 2 (cr\R(l)\ 2 ) 



¥0 



r 



Since all the terms in the previous expression but the The results derived are valid for any Bravais lattice, 
first one are analytical functions of the momentum, the Again, only the triangular lattice is of interest in this 



response function at large distances goes like 



n a /?(p) = 



2 2 -*7iT(l - 7 /2) p a p p 



(A19) 



paper. The tensor A a p^ v for the triangular lattice is |6£ 



= -74r^(^) 7/2 " 1 E G «( 1 ) G M 1 )^-7/2(^^ 2 (i)) 



(A20) 



2r( 7 /2) 



1 (- n 

4ir( 7 /2)^c J ^ 



(G Q (1)G„(1) V + G^(1)G V (1)V + G Q (1)G„(1)<W + G^(I)Gp(l)5 atf ) <^_ 7/2 (— i? 2 (l)) 

Ac 



5f7 i 7^(7 L ) 7/2 " 2 E G «( 1 ) G ^(l)G M (l)G,(l) ^^(JLjpfl)) 
81 (7/2) A c ^ A c 

_J_ ( ^ ) 7/2+2^ iiQ(1) ^ (1) ^ (1)i?i/( l ) ^ 1+7/2( JL^2 (1)) 



1/0 



c 



r 



The form of the elastic tensor can be parameterized by 
two coefficients 0(7) and 77(7) 



a result that can just follows from the symmetry proper- 
ties of the triangular lattice 62] 



^7/2 



[SpvSaP + p{.l){$Ha&vP + tinptiva] , (A21) 
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APPENDIX B: 



MATHEMATICAL IDENTITIES 
USED 



In this section useful mathematical identities are listed 
without further remarks to make the paper as self- 
contained as possible and for the purpose of fixing the 
notation. 



1. Identities present in Ewald Sums 



The Gamma identity 



1 



1 



• Misra function definition 



dtr^e-^ 2 



oo 



dtt n e - zt 



(Bl) 



(B2) 



• Misra function derivatives 
Vg</? ra (a|x — m| ) = — 2a(x — rn)Lpi +n (a\x — m| ) (B3) 



d 2 



d a d, 



-ip n (a\x— m| 2 ) = -2a5 a pip 1+n (a\x - rh\ 2 ) 
+4a 2 (f - m) a {x - m)i3(pn+2(a\x — rh\ 2 ) (B4) 



• Poisson summation formula for Gaussian integrals 



,iq-(x-R(n)) -ta\x-R(\)\ 2 



A C t(T 



,l(G+q)-x, 



k7+i-T 



(B5) 
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